import panel
import ipywidgets
panel.extension('ipywidgets')
pkg_ver = lambda pkg: "{:<20}{:}".format(pkg.__name__,pkg.__version__)
# ROOT
import uproot
print(pkg_ver(uproot))
import ROOT
# Machine Learning
import sklearn
print(pkg_ver(sklearn))
import torch
print(pkg_ver(torch))
# Data science
import scipy
print(pkg_ver(scipy))
import numpy
print(pkg_ver(numpy))
import pandas
print(pkg_ver(pandas))
import awkward
print(pkg_ver(awkward))
# Visualizations
import matplotlib
print(pkg_ver(matplotlib))
import matplotlib.pyplot as plt
import hvplot
print(pkg_ver(hvplot))
import hvplot.pandas
pandas.options.plotting.backend = 'holoviews'
import tqdm
print(pkg_ver(tqdm))
from tqdm.auto import tqdm
tqdm.pandas()
import glob
uproot 5.0.7 Welcome to JupyROOT 6.28/02 sklearn 1.2.2 torch 2.0.0 scipy 1.10.1 numpy 1.23.5 pandas 1.5.3 awkward 2.2.0 matplotlib 3.7.1 hvplot None
tqdm 4.62.3
%jsroot
"""
file = uproot.open(Data["Na22-40"](100))
c = file['COINCIDENCES;1'].arrays(library="ak")
# Gets ttree as dataframe
df = awkward.to_dataframe(c)
# Only multiplicity 2
df = df.query("detector_multiplicity == 2 & deposited_energy != 0")
# Only if it includes the scatterer of A
df = df.loc[df[df.index.get_level_values('subentry').isin([0])].index.get_level_values('entry')]
# Only if it only includes crystals of A
df = df[df.index.get_level_values('subentry').isin([0,1,2,3,4])]
# Only if 2 entries according to the previous filters
df = df.loc[df.index.get_level_values('entry')[df.index.get_level_values('entry').duplicated()]]
# Apply energy calibration to individual entries
df["deposited_energy"] = df.progress_apply(lambda x: numpy.polyval(calibrations_df.query(f"crystal == 'A{x.name[1]}' & cw == '100'").iloc[0][[0,1,2]][::-1],x.deposited_energy), axis=1)
# Recalculate sum of energies
new_sum = df.groupby(by="entry", axis=0).deposited_energy.sum()
df["total_deposited_energy"] = df.progress_apply(lambda x: new_sum.loc[x.name[0]], axis=1)
# Show dataframe
df
hist = ROOT.TH1F("","Coincidences;Energy [keV];Counts",2000,0,2000)
for i in df.iterrows():
hist.Fill(i[1].total_deposited_energy)
TH1D_draw(hist).Draw()
"""
'\nfile = uproot.open(Data["Na22-40"](100))\n\nc = file[\'COINCIDENCES;1\'].arrays(library="ak")\n\n# Gets ttree as dataframe\ndf = awkward.to_dataframe(c)\n# Only multiplicity 2\ndf = df.query("detector_multiplicity == 2 & deposited_energy != 0")\n# Only if it includes the scatterer of A\ndf = df.loc[df[df.index.get_level_values(\'subentry\').isin([0])].index.get_level_values(\'entry\')]\n# Only if it only includes crystals of A\ndf = df[df.index.get_level_values(\'subentry\').isin([0,1,2,3,4])]\n# Only if 2 entries according to the previous filters\ndf = df.loc[df.index.get_level_values(\'entry\')[df.index.get_level_values(\'entry\').duplicated()]]\n# Apply energy calibration to individual entries\ndf["deposited_energy"] = df.progress_apply(lambda x: numpy.polyval(calibrations_df.query(f"crystal == \'A{x.name[1]}\' & cw == \'100\'").iloc[0][[0,1,2]][::-1],x.deposited_energy), axis=1)\n# Recalculate sum of energies\nnew_sum = df.groupby(by="entry", axis=0).deposited_energy.sum()\ndf["total_deposited_energy"] = df.progress_apply(lambda x: new_sum.loc[x.name[0]], axis=1)\n# Show dataframe\ndf\n\nhist = ROOT.TH1F("","Coincidences;Energy [keV];Counts",2000,0,2000)\nfor i in df.iterrows():\n hist.Fill(i[1].total_deposited_energy)\n \nTH1D_draw(hist).Draw()\n'
class spectrum:
def __init__(self, File_, Configuration_, Window_, Time_):
self.__FileName = File_
self.__Configuration = Configuration_
self.__Window = Window_
self.__Time = Time_
self.__File = uproot.open(File_)
self.__DF = self.DF()
self.__TH1D = self.TH1D()
def __call__(self, ch):
return numpy.polyval(self.__Calibration[::-1],ch)
def File(self):
return self.__File
def FileName(self):
return self.__FileName
def DF(self):
# Gets ak dataframe
df_ak = self.__File['COINCIDENCES;1'].arrays(library="ak")
#Transform from ak to pandas
df = awkward.to_dataframe(df_ak)
# Only multiplicity 2
df = df.query("detector_multiplicity == 2 & deposited_energy != 0")
# Only if it includes the scatterer of A
df = df.loc[df[df.index.get_level_values('subentry').isin([0])].index.get_level_values('entry')]
# Only if it only includes crystals of A
df = df[df.index.get_level_values('subentry').isin([0,1,2,3,4])]
# Only if 2 entries according to the previous filters
df = df.loc[df.index.get_level_values('entry')[df.index.get_level_values('entry').duplicated()]]
# Apply energy calibration to individual entries
df["deposited_energy"] = df.progress_apply(lambda x: numpy.polyval(
calibrations_df.query(f"crystal == 'A{x.name[1]}' & cw == '100'").iloc[0][[0,1,2]][::-1],x.deposited_energy
), axis=1)
# Recalculate sum of energies
new_sum = df.groupby(by="entry", axis=0).deposited_energy.sum()
df["total_deposited_energy"] = df.progress_apply(lambda x: new_sum.loc[x.name[0]], axis=1)
return df
def TH1D(self):
# Create histogram
hist_s = ROOT.TH1F("","Coincidences (Sum);Energy [keV];Counts",2000,0,2000)
hist_1 = ROOT.TH1F("","Coincidences (Absorber 1);Energy [keV];Counts",2000,0,2000)
hist_2 = ROOT.TH1F("","Coincidences (Absorber 2);Energy [keV];Counts",2000,0,2000)
hist_3 = ROOT.TH1F("","Coincidences (Absorber 3);Energy [keV];Counts",2000,0,2000)
hist_4 = ROOT.TH1F("","Coincidences (Absorber 4);Energy [keV];Counts",2000,0,2000)
# Add entries
for i in self.__DF.iterrows():
match i[0][1]:
case 1:
hist_1.Fill(i[1].total_deposited_energy)
case 2:
hist_2.Fill(i[1].total_deposited_energy)
case 3:
hist_3.Fill(i[1].total_deposited_energy)
case 4:
hist_4.Fill(i[1].total_deposited_energy)
case _:
continue
hist_s.Fill(i[1].total_deposited_energy)
# Scale
hist_1.Scale(1/self.__Time)
hist_2.Scale(1/self.__Time)
hist_3.Scale(1/self.__Time)
hist_4.Scale(1/self.__Time)
hist_s.Scale(1/self.__Time)
return [hist_s, hist_1, hist_2, hist_3, hist_4]
def get_TH1D(self, i=0):
return self.__TH1D[i]
def Configuration(self):
return self.__Configuration
def Window(self):
return self.__Window
def Rate(self):
return self.get_TH1D().Integral()
calp = glob.glob('/run/media/bgameiro/d043b5e4-57a4-457e-8839-cb1adc9c72bc/Data/Calibration/**/*.CALp', recursive=True)
calibrations = []
for i in calp:
line = pandas.read_csv(
i,
sep = " : ",
skiprows=[0,4,5,6],
header=None,
engine="python"
).drop([0], axis=1).T
line["crystal"] = i.split("/")[-1][:2]
line["cw"] = i.split("/")[-1][5:8]
calibrations.append(line)
calibrations_df = pandas.concat(calibrations)
Data = {
#"Na22-0": lambda win: f"/run/media/bgameiro/d043b5e4-57a4-457e-8839-cb1adc9c72bc/Data/data_2023_05_10/Na22_center_iTEDA_moving_D.2023_05_10_T.18_10_18_C.itedABCD_lab_custom_2023.02.22_4.0v_887_900s_x0.0_y0.0_CW{win}.root",
#"Na22-10": lambda win: f"/run/media/bgameiro/d043b5e4-57a4-457e-8839-cb1adc9c72bc/Data/data_2023_05_10/Na22_center_iTEDA_moving_D.2023_05_10_T.18_25_35_C.itedABCD_lab_custom_2023.02.22_4.0v_887_900s_x10.0_y0.0_CW{win}.root",
#"Na22-20": lambda win: f"/run/media/bgameiro/d043b5e4-57a4-457e-8839-cb1adc9c72bc/Data/data_2023_05_10/Na22_center_iTEDA_moving_D.2023_05_10_T.18_40_53_C.itedABCD_lab_custom_2023.02.22_4.0v_887_900s_x20.0_y0.0_CW{win}.root",
#"Na22-30": lambda win: f"/run/media/bgameiro/d043b5e4-57a4-457e-8839-cb1adc9c72bc/Data/data_2023_05_10/Na22_center_iTEDA_moving_D.2023_05_10_T.18_56_16_C.itedABCD_lab_custom_2023.02.22_4.0v_887_900s_x30.0_y0.0_CW{win}.root",
#"Na22-40": lambda win: f"/run/media/bgameiro/d043b5e4-57a4-457e-8839-cb1adc9c72bc/Data/data_2023_05_10/Na22_center_iTEDA_moving_D.2023_05_10_T.19_11_36_C.itedABCD_lab_custom_2023.02.22_4.0v_887_900s_x40.0_y0.0_CW{win}.root",
#"Cs137-0": lambda win: f"/run/media/bgameiro/d043b5e4-57a4-457e-8839-cb1adc9c72bc/Data/data_2023_05_11/Cs137_center_iTEDA_moving_D.2023_05_11_T.09_51_47_C.itedABCD_lab_custom_2023.02.22_4.0v_887_900s_x0.0_y0.0_CW{win}.root",
#"Cs137-1": lambda win: f"/run/media/bgameiro/d043b5e4-57a4-457e-8839-cb1adc9c72bc/Data/data_2023_05_17/Cs137_mid_iTEDA_far_D.2023_05_17_T.20_57_08_C.itedABCD_lab_custom_2023.02.22_4.0v_887_1200s_serie2_2_CW{win}.root",
"Cs137-0": lambda win: f"/run/media/bgameiro/d043b5e4-57a4-457e-8839-cb1adc9c72bc/Data/data_2023_05_18/Cs137_center_iTEDA_far_D.2023_05_18_T.18_31_25_C.itedABCD_lab_custom_2023.02.22_4.0v_887_1200s_CW{win}.root",
"Cs137-10": lambda win: f"/run/media/bgameiro/d043b5e4-57a4-457e-8839-cb1adc9c72bc/Data/data_2023_05_11/Cs137_center_iTEDA_moving_D.2023_05_11_T.10_07_01_C.itedABCD_lab_custom_2023.02.22_4.0v_887_900s_x10.0_y0.0_CW{win}.root",
"Cs137-20": lambda win: f"/run/media/bgameiro/d043b5e4-57a4-457e-8839-cb1adc9c72bc/Data/data_2023_05_11/Cs137_center_iTEDA_moving_D.2023_05_11_T.10_22_17_C.itedABCD_lab_custom_2023.02.22_4.0v_887_900s_x20.0_y0.0_CW{win}.root",
"Cs137-30": lambda win: f"/run/media/bgameiro/d043b5e4-57a4-457e-8839-cb1adc9c72bc/Data/data_2023_05_11/Cs137_center_iTEDA_moving_D.2023_05_11_T.10_37_33_C.itedABCD_lab_custom_2023.02.22_4.0v_887_900s_x30.0_y0.0_CW{win}.root",
"Cs137-40": lambda win: f"/run/media/bgameiro/d043b5e4-57a4-457e-8839-cb1adc9c72bc/Data/data_2023_05_11/Cs137_center_iTEDA_moving_D.2023_05_11_T.10_52_49_C.itedABCD_lab_custom_2023.02.22_4.0v_887_900s_x40.0_y0.0_CW{win}.root",
}
def get_resolution(spectr,en):
TH1D = spectr.get_TH1D()
TH1D.GetXaxis().SetRange(TH1D.FindBin(en-250),TH1D.FindBin(en+200))
MaxBin = TH1D.FindBin(TH1D.GetMaximumBin())
ADC_Low = MaxBin-150
ADC_High = MaxBin+150
gaussFit = ROOT.TF1("gaussFit", "gaus(0)+pol3(3)", ADC_Low, ADC_High)
gaussFit.SetParameters(TH1D.GetMaximum(),MaxBin,10,0.375,-0.0005,0,0)
TH1D.Fit(gaussFit,"QR")
sigma = abs(gaussFit.GetParameter(2))
centroid = gaussFit.GetParameter(1)
return sigma*numpy.sqrt(2*numpy.log(2))*2/centroid*100, centroid, abs(centroid-en)/en*100
def TH1D_draw(spectr,en):
canvas = ROOT.TCanvas()
canvas.cd()
TH1D = spectr.get_TH1D()
#TH1D.SetTitle(repr(cell))
TH1D.SetStats(False)
latex = ROOT.TLatex()
latex.SetNDC()
latex.SetTextSize(0.03)
TH1D.Draw("PE PLC")
l1,l2,l3 = get_resolution(spectr,en)
latex.DrawText(0.7, 0.8, "R: {:.2f}%".format(l1))
latex.DrawText(0.7, 0.7, "E: {:.0f}keV".format(l2))
TH1D1 = spectr.get_TH1D(1)
TH1D1.Draw("SAME PLC")
TH1D2 = spectr.get_TH1D(2)
TH1D2.Draw("SAME PLC")
TH1D3 = spectr.get_TH1D(3)
TH1D3.Draw("SAME PLC")
TH1D4 = spectr.get_TH1D(4)
TH1D4.Draw("SAME PLC")
legend = ROOT.TLegend(0.1,0.7,0.9,0.9)
legend.AddEntry(TH1D,"Sum of absorbers","l")
legend.AddEntry(TH1D1,"Absorber 1","l")
legend.AddEntry(TH1D2,"Absorber 2","l")
legend.AddEntry(TH1D3,"Absorber 3","l")
legend.AddEntry(TH1D4,"Absorber 4","l")
legend.Draw("SAME")
canvas.BuildLegend()
return canvas
entries = []
for i in tqdm(Data):
source, conf = i.split("-")
for CW in tqdm([100,150,200,250]):
spectr = spectrum(
Data[i](CW),
int(conf),
CW,
900 if "900s" in Data[i](CW) else 1200
)
en = 662 if source == "Cs137" else 511
res, cent, fit = get_resolution(spectr,en)
entries.append(
pandas.DataFrame({
"resolution": res,
"fit": fit,
"cps": spectr.Rate(),
"cw": CW,
"configuration": int(conf),
"obj": spectr,
}, index=[0])
)
entries_df = pandas.concat(entries, ignore_index=True)
0%| | 0/5 [00:00<?, ?it/s]
0%| | 0/4 [00:00<?, ?it/s]
0%| | 0/68550 [00:00<?, ?it/s]
0%| | 0/68550 [00:00<?, ?it/s]
0%| | 0/61970 [00:00<?, ?it/s]
0%| | 0/61970 [00:00<?, ?it/s]
0%| | 0/58044 [00:00<?, ?it/s]
0%| | 0/58044 [00:00<?, ?it/s]
0%| | 0/55510 [00:00<?, ?it/s]
0%| | 0/55510 [00:00<?, ?it/s]
0%| | 0/4 [00:00<?, ?it/s]
0%| | 0/52188 [00:00<?, ?it/s]
0%| | 0/52188 [00:00<?, ?it/s]
0%| | 0/46988 [00:00<?, ?it/s]
0%| | 0/46988 [00:00<?, ?it/s]
0%| | 0/43842 [00:00<?, ?it/s]
0%| | 0/43842 [00:00<?, ?it/s]
0%| | 0/42110 [00:00<?, ?it/s]
0%| | 0/42110 [00:00<?, ?it/s]
0%| | 0/4 [00:00<?, ?it/s]
0%| | 0/77886 [00:00<?, ?it/s]
0%| | 0/77886 [00:00<?, ?it/s]
0%| | 0/69772 [00:00<?, ?it/s]
0%| | 0/69772 [00:00<?, ?it/s]
0%| | 0/64626 [00:00<?, ?it/s]
0%| | 0/64626 [00:00<?, ?it/s]
0%| | 0/61522 [00:00<?, ?it/s]
0%| | 0/61522 [00:00<?, ?it/s]
0%| | 0/4 [00:00<?, ?it/s]
0%| | 0/97626 [00:00<?, ?it/s]
0%| | 0/97626 [00:00<?, ?it/s]
0%| | 0/87606 [00:00<?, ?it/s]
0%| | 0/87606 [00:00<?, ?it/s]
0%| | 0/81152 [00:00<?, ?it/s]
0%| | 0/81152 [00:00<?, ?it/s]
0%| | 0/77240 [00:00<?, ?it/s]
0%| | 0/77240 [00:00<?, ?it/s]
0%| | 0/4 [00:00<?, ?it/s]
0%| | 0/121068 [00:00<?, ?it/s]
0%| | 0/121068 [00:00<?, ?it/s]
0%| | 0/108498 [00:00<?, ?it/s]
0%| | 0/108498 [00:00<?, ?it/s]
0%| | 0/100654 [00:00<?, ?it/s]
0%| | 0/100654 [00:00<?, ?it/s]
0%| | 0/95422 [00:00<?, ?it/s]
0%| | 0/95422 [00:00<?, ?it/s]
Info in <TCanvas::MakeDefCanvas>: created default TCanvas with name c1
entries_df["resolution"] = entries_df.obj.map(lambda x: get_resolution(x,662)[0])
entries_df["fit"] = entries_df.obj.map(lambda x: get_resolution(x,662)[2])
entries_df
| resolution | fit | cps | cw | configuration | obj | |
|---|---|---|---|---|---|---|
| 0 | 8.940606 | 2.044403 | 19.550833 | 100 | 0 | <__main__.spectrum object at 0x7f434d1fb820> |
| 1 | 8.883976 | 1.959143 | 17.604167 | 150 | 0 | <__main__.spectrum object at 0x7f434a2d3a60> |
| 2 | 8.816424 | 1.898117 | 16.409167 | 200 | 0 | <__main__.spectrum object at 0x7f434a275a20> |
| 3 | 9.137346 | 1.949307 | 15.675833 | 250 | 0 | <__main__.spectrum object at 0x7f434878a950> |
| 4 | 8.636753 | 3.762054 | 19.307778 | 100 | 10 | <__main__.spectrum object at 0x7f434a3eb520> |
| 5 | 8.574016 | 3.721391 | 17.285556 | 150 | 10 | <__main__.spectrum object at 0x7f43482b1c90> |
| 6 | 8.491103 | 3.632212 | 16.070000 | 200 | 10 | <__main__.spectrum object at 0x7f42d1f7ae30> |
| 7 | 8.463257 | 3.594018 | 15.391111 | 250 | 10 | <__main__.spectrum object at 0x7f434a3eb640> |
| 8 | 8.926223 | 4.546599 | 29.617778 | 100 | 20 | <__main__.spectrum object at 0x7f42d13ed5d0> |
| 9 | 8.420145 | 4.330549 | 26.447778 | 150 | 20 | <__main__.spectrum object at 0x7f433bfd8250> |
| 10 | 8.378418 | 4.257623 | 24.390000 | 200 | 20 | <__main__.spectrum object at 0x7f4349818d90> |
| 11 | 8.314041 | 4.196396 | 23.146667 | 250 | 20 | <__main__.spectrum object at 0x7f42d08597b0> |
| 12 | 8.873971 | 5.154211 | 37.441111 | 100 | 30 | <__main__.spectrum object at 0x7f42c9ee2b30> |
| 13 | 8.794574 | 5.072420 | 33.440000 | 150 | 30 | <__main__.spectrum object at 0x7f42d1972a70> |
| 14 | 8.721396 | 4.970308 | 30.840000 | 200 | 30 | <__main__.spectrum object at 0x7f42cf450820> |
| 15 | 8.715796 | 4.923926 | 29.320000 | 250 | 30 | <__main__.spectrum object at 0x7f42ce866350> |
| 16 | 8.833005 | 4.799102 | 47.157778 | 100 | 40 | <__main__.spectrum object at 0x7f42d0c39270> |
| 17 | 8.804763 | 4.706591 | 42.051111 | 150 | 40 | <__main__.spectrum object at 0x7f42d02ed4e0> |
| 18 | 8.805299 | 4.643221 | 38.850000 | 200 | 40 | <__main__.spectrum object at 0x7f434a17fc70> |
| 19 | 8.819354 | 4.616517 | 36.714444 | 250 | 40 | <__main__.spectrum object at 0x7f42d18722f0> |
entries_df.hvplot.scatter("configuration","resolution",by="cw")
entries_df.hvplot.scatter("configuration","fit",by="cw")
entries_df.hvplot.scatter("configuration","cps",by="cw")
l=entries_df.loc[0].obj
get_resolution(l,662)
TH1D_draw(l,662).Draw()
l=entries_df.loc[16].obj
get_resolution(l,662)
TH1D_draw(l,662).Draw()